Diffusion and jump-length distribution in liquid and amorphous Cu33Zr67 
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Using molecular dynamics simulation, we calculate the distribution of atomic jumps in Cu33Zre7 
in the liquid and glassy states. In both states the distribution of jump lengths can be described 
by a temperature independent exponential of the length and an effective activation energy plus a 
contribution of elastic displacements at short distances. Upon cooling the contribution of shorter 
jumps dominates. No indication of an enhanced probability to jump over a nearest neighbor distance 
was found. We find a smooth transition from flow in the liquid to jumps in the glass. The correlation 
factor of the diffusion constant decreases with decreasing temperature, causing a drop of diffusion 
below the Arrhenius value, despite an apparent Arrhenius law for the jump probability. 

PACS numbers: 61.43.Fs,66.10.-x,66.30.-h,64.70.Pf 
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I. INTRODUCTION 

Metallic glasses are of great interest both fundamen- 
tally and due to numerous applications. The absence 
of strong covalent bonds makes them the prime example 
of random dense packing. Atomic transport in metallic 
glasses and their melts is generally considered to be ef- 
fected by two distinct mechanisms: flow in the liquid, 
and thermally activated hopping in the glassy state. See 
Ref. [l| for a review. When a glass forming liquid is 
cooled towards the glass transition diffusion drops faster 
than predicted by an Arrhenius law. The diffusion coef- 
ficient is often described by a Vogel-Fulcher law which 
can be rationalized e. g. by a free volume descrip- 
tion. Mode coupling theory predicts for fragile glass 
formers, to which amorphous metallic alloys belong, a 
dynamical transition at a critical temperature, T c , well 
above the glass transition temperature T g . At T c flow 
motion freezes and diffusion vanishes with a power law 
D on (T — Tc) 1 apart from a residual contribution from 
atomic hopping-. The mechanism of this change from 
flow to hopping is not understood. 

In liquids with loose packing, flow is governed by bi- 
nary collisions between the constituent atoms. In the un- 
dercooled melt, where the atomic packing is dense, flow 
is strongly collective. This is reflected in the observed 
small isotope effect of diffusion^, similar to the one ob- 
served earlier in the glassy stated. From these experi- 
ments it has been concluded that, both below and above 
the glass transition, diffusion is by collective motion of 
ten or more atoms. This poses the question of whether 
the change from flow to hopping is a change to a new el- 
ementary process or whether hopping evolves out of the 
flow motion. 

Over the last few years computer simulations have pro- 
vided considerable insight into the atomic dynamics of 
glasses and undercooled liquids. Early molecular dy- 
namics simulations (MD) have shown collective jumps in 
undercooled liquids^. Chains of atoms replacing each 
other were observed, i.e. the single atoms jumped by a 
nearest neighbor distance. These chains can close to form 
rings^i Comparing successive configurations, averaged 



over typical vibrational times, one again finds chain-like 
structures of atoms which have moved collectively in the 
undercooled liquid 8-910 . These are not necessarily re- 
placement chains. The jump process in the glassy state 
has been found to involve many atoms, each single atom 
moving only a fraction of the nearest neighbor distance 
in such a jumpAi*i2iii 

Upon cooling towards the glass-transition a striking 
feature is seen in the self-part of the van Hove function 
G s (r, t) which is related to the probability that an atom 
has moved by a distance r during a time t. At high 
temperatures G s (r, t) is perfectly Gaussian and broad- 
ens oc \fi. Upon cooling towards T c , and beyond, a tail 
to larger distances grows with time. Finally approach- 
ing T c an additional second peak at the nearest neigh- 
bor distance evolves, particularly for the more mobile 
components. This effect can be taken as one of the sig- 
natures of the glass transition^. From this behavior it 
was concluded that there is a single peaked distribution 
of hopping distances^. The time evolution of G s (r, i) 
in CuZr could be reproduced by a simple model involv- 
ing jumps over nearest neighbor distances plus a residual 
small flow 16 . In this picture the jump motion dominates 
the diffusion in the undercooled liquid and the super- 
Arrhcnius drop of diffusion stems from an increase of the 
return jump probability, as one would expect from an 
increasing number of blocked paths. 

In a quantitative investigation of the deviation of 
G s (r,t) from a Gaussian, the non-Gaussianity parame- 
ter, a,<2,{T,t), was found to increase rapidly in the under- 
cooled liquid but no abrupt change near T c was seeniiA£. 
The time evolution could be understood from a model of 
collective jumpsit 

Inspecting the pressure derivative of the diffusion con- 
stant, the apparent activation volume, one finds a strong 
cusp at T c which could indicate a change of diffusion 
mechanismiS. 

These different findings pose the question of whether 
there is a change in the elementary process of diffusion 
near T c , in particular whether one might observe the evo- 
lution of a typical jump process in the glass. 
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II. SIMULATION DETAILS 

Here we report a MD investigation of the atomic jump 
lengths in undercooled and glassy Zr 6 7Cu33. For the 
inter-atomic interaction we use a modified Embedded 
Atom Method. 20 The parameters were fitted to repro- 
duce the experimental values of Cu, Zr and CuZr2 crys- 
tals. The universal energy-volume relation of Rose et 
alm^ was used to determine the anharmonic contribu- 
tions, not sampled in the crystal but of essential im- 
portance in the disordered glassy state. We get lat- 
tice parameters a — 0.363, a = 0.323, /, c = 0.516 
and a — 0.338, c = 10.35 nm (experimental values^ 
a = 0.362, a = 0.323, c = 0.515 and a = 0.322, 
c = 11.18 nm) for Cu, Zr and CuZr2, respectively. The 
CuZr 2 lattice is slightly distorted. The atomic volume, 
however, is only 2% too large. The sublimation ener- 
gies for Cu and Zr (3.53 and 6.34 eV) agree with exper- 
iment. We find enthalpies of fusion per atom relative to 
the mono-atomic crystalline phases at room temperature 
of A f H = 0.22 and 0.18 eV (experiment for CuZr 2 and 
CuZr, respectively. The vacancy formation energies are 
1.32 and 1.63 eV (experiment 1.28 and > 1.5 eV) for 
Cu and Zr, respectively. Additionally the phonon dis- 
persion curves and elastic constants of the mono-atomic 
lattices were used. In the case of Cu excellent agreement 
was achieved. In Zr we get an overall agreement with 
experiment but some phonons deviate up to 30%, sim- 
ilar to other worltt. No attempt was made to fine 
tune the potential to fully reproduce the phase diagram. 
The detailed form and the parameters are given in the 
appendix. For more details on the fitting procedure see 
Reft. 

The MD calculations were done using the velocity Ver- 
let algorithm with a time-step of 2.5 TO -15 s and systems 
of N = N Zl + N Cu = 1000 atoms with periodic bound- 
ary conditions. Previous work on other systems (e. g. 
soft spheres, binary Lennard-Jones, Se) by us and other 
groups has shown that this size suffices to reproduce the 
dynamics at elevated temperatures. As additional test 
some runs with N — 8000 were done for comparison. For 
the questions investigated in this work long aging times 
are more important than large system sizes. The pres- 
sure was kept constant following Ref. 28 ^ using a volume 
mass of w \/N ■ mz r and an additional damping term 
to prevent oscillations. Temperature was controlled by a 
Nose-thermostat following Hoover—. 

Three independent samples were prepared by a quench 
from the hot liquid and were aged in intermittent stages, 
as shown in Fig. ^ We cool, in steps of 100 K, with a 
rate of 10 12 K/s from 2000 K to the simulation tempera- 
ture. At each temperature step the samples were aged for 
times ranging from 1 ns at 2000 K to 2 ns for T < 1000 K, 
before continuing the quench. The effective cooling rate 
was thus lowered by about an order of magnitude, com- 
pared to a straight quench with a constant rate. Before 
the actual measurements at a given temperature the sys- 
tems were additionally aged at constant T for different 



times, up to 5.5 ns, as indicated by the dotted lines in 

Fig.m 
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FIG. 1: Quench history of the samples (solid line) . The dotted 
lines show the lengths of additional aging before the start of 
the measuring runs. 



III. GLASS TRANSITION AND AGING 

Experimentally Zrg7Cu33 is a good glass-former, which 
melts at 1310 K 30 and can be under-cooled by a few hun- 
dred Kelvin. The experimental glass-transition temper- 
ature varies from 600 K at cooling-rates of 0.01 K/s to 
750 K at 10 6 K/t. To determine the glass transition 
temperature of our CuZr model we did additional runs 
without the intermittent aging and monitored the poten- 
tial energy and the volume as function of temperature for 
different quench rates, ranging from Q = 2.5 • 10 13 K/s 
to Q = 4 • 10 10 K/s, see Fig. |2 for the potential energy. 
The glass transition temperature, T g , was defined by the 
crossover from the low temperature to the liquid behav- 
ior. In the limited range of Q covered by the simulation 
the dependence of T g on Q can be expressed by a loga- 
rithmic law^i 



T g = (463 + 20.6- In Q) K. 



(1) 



For our lowest quench rate Q = 4 • 10 10 K/s we find 
T g = 965 K extrapolating to Q = 4 • 10 6 K/s Eq.UJgives 
Tg = 747 K in excellent agreement with the experimental 
value. 

From Fig.[2]together with Eq.^we nud a linear depen- 
dence of the average potential energy per atom at T = 
on the glass transition temperature 

E P {T = 0, Q) = (-5.7715 + 1.5 • 10~ 4 T g (Q)) eV. (2) 

Such a linear equation was also observed in amorphous 
Se for both energy and atomic volume 32 . In CuZr the 
quench rate dependence of the volume is too small to be 
evaluated. 
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FIG. 2: Average potential energy per atom as function of 
temperature for different quench rates: from top to bottom 
Q = 2.5 • 10 13 , 5 ■ 10 12 , 1 ■ 10 12 , 2 ■ 10 11 and 4 • 10 10 K/s. 



Monitoring the average potential energy per atom at 
T = 1000 K, near the glass transition temperature, start- 
ing at the end of the quench, Fig. ^ we find that it drops 
in 35 ns by about 0.014 ev. The statistics is not sufficient 
to determine the decay law unambiguously. 




10 10 
time (s) 



FIG. 3: Average mean square displacement as function of 
time on a double logarithmic scale. Temperature from top to 
bottom: 2000, 1800, 1600, 1400, 1200, 1100, 1000, 960, 900, 
800, 700 K. 



IV. DIFFUSION AND HETEROGENEITY 



The average atomic mean square displacements 
(msqd), from the respective configurations after aging, 
are shown in Fig.[31in a double logarithmic plot. For short 
times one observes an increase oc t 2 which is typical for 
vibrational and ballistic motion. For long times the msqd 
increases oc t, indicative of long range diffusion. Lower- 
ing the temperatures below T = 1400 K the plateau, 
typical for the undercooled liquid and the glass, evolves 
between these two limits. This onset of the plateau cor- 
rellates well with the experimental melting temperature 
T m = 1310 KS1. At the lowest temperatures, in the ps 
range, one can just see some small wiggles which reflect 
the vibration spectrum. 

The diffusion coefficients of the two components, 
Fig. 01 were calculated in the usual way from the slope 
of the long time limit of the msqd. Fitting the diffusion 
coefficients in the undercooled melt according to MCT, 



£>mct(T) = j D mct /(T- 



(3) 



we find T c = 1025 K and 7 = 1.92 and 1.34 for Zr and 
Cu, respectively. A fit to the relaxation time of the in- 
termediate scattering function with the same T c gives 
slightly higher values 7 = 2.2 and 1.57 for Zr and Cu, 
respectively 1,33 . These numbers are meant as a guide to 
the relevant temperatures and are not exact. An alterna- 
tive fit with the Volgel-Fulcher-Tammann (VFT) relation 

D VFT (T) = D^ T exp (~E VFT /k(T - T VFT )) (4) 
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FIG. 4: Diffusion coefficients in Zr67Cu33 (Zr: diamonds, Cu: 
spheres) The dashed lines represent a fit with MCT, using the 
same temperature T c for both components, in the undercooled 
melt and an Arrhenius fit in the glass. 



gives for Zr D^ FT = 2.65 • 10" 8 m 2 /s, E VFT = 0.28 eV, 
T VFT = 641 K and for Cu L# FT = 3.70 • 10~ 8 m 2 /s, 
£ VFT = 0.30 eV, T VFT = 477 K. In the glass the diffu- 
sivity can be described by the usual Arrhenius law 

D A rrh(T) = Dfr" h exp (-E m /kT). (5) 

Neglecting the values at T = 800 K, which are proba- 
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bly too high due to too short aging, we obtain Dg = 
1.09 • 1CT 3 m 2 /s, E m = 1.7 eV, for Zr and D$ Ilh = 
1.41 • 1CT 5 m 2 /s, E m = 1.1 eV, for Cu. Due to the small 
fitting range there is a considerable margin of error on 
these values. The deviation of the present values from the 
ones reported earlier by GaukelSI is due to the much im- 
proved statistics and aging of the present work. It should 
further be noted that fits in the undercooled liquid using 
the VFT or MCT expression is strongly influenced by 
the range of temperatures included in the fit. Further- 
more in the fits of the previous work different values of 
T c were allowed for Cu and Zr, respectively, whereas in 
the present work the condition of a unique value T c was 
imposed. The data do not suffice to validate or invalidate 
this condition. 

Including the 800 K values the activation energies 
would be considerably smaller (E m = 0.99 and 0.75 eV 
for Zr and Cu, respectivelji^. 

Our results agree well with simulations of the simi- 
lar NiZr system where a totally different model for the 
inter-atomic interaction was used— The diffusion coeffi- 
cients of the two components both in NiZr and in a binary 
Lennard- Jonesi^ glass at zero pressure are nearly parallel 
in the melt, whereas in Zr67Cu33 they diverge. This effect 
is probably due to the weaker coupling between the two 
components in CuZr. Experimentally this is reflected in 
lower enthalpy of fusion in CuZr 23 compared to NiZr—. 

In isotropic diffusion the atomic displacements are 
Gaussian distributed. In undercooled liquids and in 
glasses Gaussianity is violated over long time scales. This 
non-Gaussianity indicates different mobilities of differ- 
ent atoms over long time scales. This so called dy- 
namic heterogeneity is quantified by the non-Gaussianity 
parameter— 

3 < Ar\t) > , . 

°*® = 5<Ar»(t) > 2 X ' (6) 



where < ... > denotes time averaging and Ar 2 (t) and 
Ar 4 (£) are the mean square and quartic displacements. 

Fig.[S]shows this non-Gaussianity in the Cu-subsystem 
for temperatures from 800 to 2000 K. The general be- 
havior resembles the one observed in other glass form- 
ers. In the liquid above 1400 K dynamic heterogeneity 
is weak and due to different local vibrational densities of 
state. The maximal non-Gaussianity is at typical vibra- 
tional times (ps) . With increasing undercooling and even 
more in the glassy state the non-Gaussianity rapidly in- 
creases. The maximum is reached later and later. Com- 
paring with Fig. |31 one sees that the decay of the non- 
Gaussianity correlates with the onset of the diffusional 
part of the msqd following the plateau. The increase of 
<x<2 (t) from its vibrational value follows the yfi law which 
was attributed to collective jump motion 17 . 
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FIG. 5: Non-Gaussianity versus time of the Cu-subsystem in 
Zr67Cu33 on a logarithmic scale. Temperatures from top to 
bottom: 800, 900, 1000, 1100, 1200, 1400, 1600, 1800 and 
2000 K. 

V. JUMP LENGTH DISTRIBUTION 

After aging we started the observation of jump pro- 
cesses. To get sufficient statistics we used observation 
times for the detection of jumps from 1.25 ns at 1400 K 
to 2.5 ns below 1200 K. As definition of jump we use a 
rapid transit of an atom between two sites of residence. 
One can see from a direct inspection of the particle tra- 
jectories that in dense liquids or glasses such "jumps" 
are. not ballistic. For the CuZr system such a trajectory 
has been shown for a jump of a Cu-atom over a near- 
est neighbor distance in Refs. Pl37j . Furthermore the 
amplitudes of short time excursions are very large which 
makes an unambiguous definition of jumps difficult. 

One method to detect atomic jumps is to study the 
"inherent dynamics" by quenching the system at given 
time intervals "instantaneously" to K and then study 
the properties of the "energy landscape" . Here we use a 
more direct approach and define jumps by differences be- 
tween atomic positions averaged over vibrational times. 
This means we define jumps of an atom in terms of abso- 
lute coordinates and not relative to the neighbors of the 
atom. 

As a first step we define for each time step and each 
atom an average atomic position by 

(R»(i)> = T / R n (t')dt'. (7) 

Since this averaging is done for each time step (R n (t)) 
is still, on the scale of the discretization by the simula- 
tion time step, a continuous function but is smoothed by 
averaging over the vibrations. An instantaneous jump 
would cause a steep ramp in (R n (t)). 

For the detection of a jump we use two criteria. First, 
the instantaneous position of an atom must differ by a 
minimum cutoff length from the average taken at a pre- 
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vious time: 



R(i)-<R*(i--*)) 



> n. 



(8) 



We want to exclude from our jump detection those excur- 
sions where an atom has a large amplitude momentarily, 
but immediately returns to its old site. Therefore, when 
the above condition is fulfilled we additionally compare 
average positions separated by the fixed time interval <5 2 

A(R n (i)) = (R"(t + 0.5<5 2 )) - (R"(i - 0.5tf 2 )>. (9) 

A jump at time to is recorded when this difference for the 
first time exceeds a limit r 2 

|A(R n (i„))| >r 2 . (10) 

The corresponding atomic jump length is defined as 

£ = A(R"(i )). (11) 

The time interval between the start time of the aver- 
aging for the final configuration, i+0.5<5 2 — 0.5<5i, and the 
end time of the averaging for the initial configuration^ — 
0.5S 2 + 0.55i is 



= 82- Si. 



(12) 



This method is rather CPU-time consuming since at 
each time step two averages have to be done for each 
atom. This is however necessary if one wants to get the 
necessary time resolution. 

Throughout the simulation we used the parameters 
Si = 2.5 ps, 82 = 4 ps and r 2 = 0.05 nm. 

The necessity of averaging the atomic configurations 
over given time intervals limits the accuracy of the 
method. If a jump is completed in the time interval, 
iwait ! its length will be measured correctly by Eq. A 
longer "jump time" leads to a reduced apparent jump 
length. On the other hand a movement of the atom with 
constant velocity v during the time <5i -I- 82 would be de- 
tected as jump if v ■ (Si + 82) > f 2 - Since averaging over 
typical vibrational times implies 8\ ~ ps we do not expect 
this to be an important limitation. 

More serious is the limited resolution of rapid succes- 
sive jumps. If two successive jumps are completed within 
the time interval t W ait they will be regarded as a single 
jump. If, however, the second jump happens near the 
end of the time interval it will be counted only partially. 
For example, should the second jump be the reverse of 
the first jump we might under some circumstances still 
record some shorter range jump, given by the sum of the 
forward jump and a fraction of the back-jump. On the 
other hand, if the first jump is rapidly followed by a suc- 
cessive forward jump, normally correlated with the first 
jump, we record an effective jump length which is too 
short. These two effects should approximately cancel for 
the investigated temperatures. For lower temperatures 
where the fraction of back-jumps increases it leads to an 



overestimate of the mean square displacement calculated 
from the jumps, compared to the exact value. 

Once the jump criteria hold for a time step they will 
normally hold for a subsequent time interval. To avoid 
double counting, therefore after each jump, we introduce 
a dead time <5 2 during which jumps are not counted. This 
will lead to small inaccuracies, mainly for the shortest 
jump lengths. 
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FIG. 6: Atomic jump rates, Tj^mp (^cut ) , on a logarithmic scale 
over temperature for Zr for different cutoff lengths, ^cut: Zr 
(top) and Cu (bottom). Lines from top to bottom: £ cut = 
0.06,0.1,0.15,0.2,0.25 nm. 



From the recorded jumps we can calculate how often 
and how far each atom jumped during the observation 
time. This defines the average atomic jump rates 
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E < rop (^,Us) (13) 



where N? u (T, £, i bs) is the number of jumps of atoms 
of species j with jump length in the interval [£ — 8£/2, £ + 
8£/2] with tabs and T the observation time and temper- 
ature, respectively. Fig. shows the rates for both com- 
ponents versus temperature for different cutoff length, 
£ < ^cut- As to be expected the jump rates increase with 
temperature. No obvious break is apparent at T c or T g . 
For the lower temperatures the jump rate for Zr is clearly 
much smaller than that for Cu, particularly for the larger 
£ cu t- This is in accordance with the lower diffusivity of 
Zr. Jumps over nearest neighbor distances are observ- 
able for Zr only well above T g . For Cu they are observed 
at all temperatures down to 800 K. At the highest tem- 
perature the rate for all jumps with £ < £ cut = 0.06 nm 
approaches the theoretical limit of resolution of this sim- 
ulation, C - t = 2.5 • 10 11 s _1 . This limits the applicabil- 
ity of the method to higher temperatures. It reflects the 
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gradual transition from jumps to flow. It becomes mean- 
ingless to assign average starting positions to the atoms. 
In Fig. El average jump rates are shown. Individual atoms 
will have higher rates. However, the dynamic heterogene- 
ity is no longer so important at T — 1400 K so that the 
deviations from the average are not too large. We will 
argue further down that the majority of apparent jumps 
with short lengths is actually due to the elastic displace- 
ment of surrounding atoms, accompanying all jumps. 

For a more detailed investigation, we recorded for each 
detected jump the time of the jump and the initial and 
final positions of the jumping atom, according to Eq. [7| 
From these data the probability that an atom jumps by a 
certain distance is calculated. Summing over the atomic 
jumps we also calculate the atomic mean square displace- 
ments which can be compared with the ones gained di- 
rectly in the simulation. 

The probability that an atom jumps over a distance I 

is 

P j (T,i)Si = ^-4 mp (T,^ obs ). (14) 

Fig.[7]shows the distribution of the numbers of jumps per 
second against jump length for the larger majority com- 
ponent Zr for temperatures ranging from 900 to 1400 K 
(top) and for the smaller minority component Cu for tem- 
peratures ranging from 800 to 1400 K (bottom). Looking 
first at the curves for Zr one clearly sees no indication of 
a preferred jump length. The distribution can be fitted 
with a simple form 

iWT, I) = A inmp e- E ^/ kT e- 1 ' 1 — (15) 

with Af* mp = 1.83 • 10 28 l/(m-s), £f u r mp = 1.51 eV and 
!? r mp = 0.033 nm. The apparent activation energy E?* mp 
agrees within some 10% with the diffusional one in the 
glassy state. In the expression for the undercooled melt, 
Eq. it corresponds to the apparent activation energy 
at 1140 K. As shown by the dotted line this fit works 
well in the whole temperature range investigated which 
spreads over both T g and T c . Of course S? r mp has to 
be interpreted as an effective activation energy. There 
will be a spread of activation energies which is absorbed 
by the prefactor ^4? u r mp . The probability of jumps over a 
nearest neighbor distance is two orders of magnitude less 
than the one for jumps over half that distance. 

In the case of Cu the situation seems more complicated. 
The probability of jumps over a nearest neighbor distance 
is only one order of magnitude less than the one for jumps 
over half that distance. Different from Zr, there is a 
distinct tail to higher jump lengths whose origin is not 
absolutely clear. One would expect a cutoff of A(R n (t)) 
for distances not much greater than the nearest neighbor 
distance. In a densely packed material a jump over larger 
distances is highly improbable. We rather think that 
the tail is due to jumps which follow so rapidly that we 
cannot resolve them. This view is supported also by the 
correlation factors exceeding unity. We will see below 




jump length (nm) 

FIG. 7: Distribution of jumps/second over jump length for 
Zr (top) and Cu (bottom). Temperatures from top to bot- 
tom: 1400, 1200, 1100, 1000 and 900 K. The dotted lines 
indicate the fits by exponential jump length distributions, see 
text. The respective nearest neighbor distances for the two 
components are indicated by the vertical dotted lines. 

that the contribution of these "extra long jumps" drops 
at T c to less than 20% of the total atomic mean square 
displacement. 

For short jump distances we see again the "elastic 
tail" which evolves into a shoulder before it merges 
with the curve for the more important long jumps. For 
these different fitting laws apply in the glass and in the 
undercooled liquid. In the glass, the Arrhenius form, 
Eq. ^3 works just as for Zr. The fitted values are: 
^Po U „ gj u mp = 9.79 • 10 26 l/(s-m), £S gjump = 0.82 eV and 
^kmgjump = 0.024 nm. This apparent activation energy is 
about 15% less than the diffusional E m used in Fig.^Jbut 
is higher than the value gained after including T — 800 K 
in the fit of the diffusivity. 

This fit breaks down in the undercooled liquid regime. 
The slopes of the curves in Fig. decrease with temper- 
ature, at variance with Eq. 1151 Extending the fit from 
the glass into the liquid one finds still quite good values 
for jumps over nearest neighbor distances. However for 
I = 0.4 nm at T = 1400 K, Ea.HKlalreadv underestimates 
the jump probability by nearly an order of magnitude. 
The jump length distribution for Cu in the undercooled 
liquid can be fitted by an exponential law with a temper- 
ature dependent ^iongjump(r) 

nCu i (T f)= BP U ■ p- £ /^longj„mp(T) J, > rj, 

longjumpV ' / longjump c 

(16) 

with £C gjump = 1.96 • 10 22 1/sm and 4>n g j um p(T) = 
0.038, 0.043 and 0.05 nm for T = 1100, 1200 and 1400 K, 
respectively. No direct physical origin for this relation is 
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evident. From the fact that the Arrhenius like tempera- 
ture dependence, Eq.^J is still applicable for jumps over 
nearest neighbor distances the most probable explanation 
seems to be that with increasing temperature more and 
more rapidly succeeding jumps, which cannot be resolved 
into separate single jumps, occur with predominantly ad- 
ditive direction. This is in accordance with the increase 
in correlation factor discussed later on in this paper. 

For both components we observe a rapid increase of 
the jump length distribution toward short jump length. 
There is a clear excess above the value given by Eq. IT51 
For Zr at T = 900 K, where the feature is most promi- 
nent, it can be fitted by Pj U m P (T,£) <x l/£ 5 , close to the 
oc 1 /^-dependence following from the 1/r decay of elas- 
tic displacements. We therefore conclude that this rapid 
increase reflects the elastic displacements accompanying 
all jumps. Due to dense packing they will also be present 
in the melt but are much more prominent in the glassy 
state. Two effects will strongly reduce their contribution 
to diffusion. First, the elastic displacement patterns are 
different for each jump leading to a cancellation of the 
bulk of these displacements after a few jumps. In the 
present system this is clearly observed in the glass when 
the Cu-atoms are by more than an order of magnitude 
more mobile than the Zr ones. The Cu-atoms move in a 
sluggish matrix of Zr-atoms. Jumps of the Cu-atoms will 
be accompanied by displacements of the Zr-atoms with- 
out destroying their topology. After a few subsequent 
jumps, dominated by the faster Cu, the Zr will more or 
less be back to their original sites. Our algorithm will 
pick up the occasional large displacements of the Zr but 
not their return in smaller steps. We will see this effect, 
further down, in the mean square displacement of Zr at 
T < 900 K, Fig. ^3 temperatures the dynamic hetero- 
geneity, meaning that at any given time only a small sub- 
set of atoms is mobile, will increase this effect and spread 
it to the second component. Secondly, after a jump pro- 
cess the glass is locally excited and relaxes toward the 
local equilibrium. This is done by comparatively slight 
shifts of the average atomic positions. This again reduces 
the contribution of the small lengths to diffusion. 

VI. AVERAGE MEAN SQUARE 
DISPLACEMENTS 

To check the importance of the different jump lengths 
for diffusion we compare the average atomic mean square 
displacement (msqd) taken directly from MD with the 
one obtained by adding the jumps used to obtain the 
distribution of Fig. [7| 

|Ai? jump (i,4ut)| 2 = ( | ]T AR"(i')| 2 )„ (17) 

*'<t 

where (. . .)„ indicates averaging over atoms and configu- 
rations. We do this for different lower cutoffs of the jump 
length, £ cut . In the undercooked melt, Figs.[S]and[!|| there 



is an excellent agreement between the exact curves and 
the ones gained this way. In the figures, we added to 

l 1 2 

I A-Rju mp (t) ^cut)| the vibrational msqd which can be ob- 
tained from the short time behavior. 
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FIG. 8: Average atomic mean square displacement versus 
time at T — 1400 K (full line). Atomic mean square dis- 
placement calculated from the jumps used for the distribution, 
Fig. 1, for different cutoffs: from top to bottom all jumps with 
jump lengths greater than 0.06, 0.1, 0.15, 0.2 and 0.25 nm. 



0.06 
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FIG. 9: Average atomic mean square displacement versus 
time at T — 1000 K (full line). Atomic mean square dis- 
placement calculated from the jumps used for the distribution, 
Fig. 1, for different cutoffs: from top to bottom all jumps with 
jump lengths greater than 0.06, 0.1, 0.15, 0.2 and 0.25 nm. 

We find for Zr that even at the temperature T — 
1400 K, i.e. about 40% above T c , the msqd and, there- 
fore also diffusion, is dominated by jumps much shorter 
than the nearest neighbor distance, R^ N ~ 0.3 nm. 
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Jumps of more than 0.25 nm contribute about 5%. At 
T = 1000 K their contribution vanishes. In Cu, jumps 
over Rpfff ~ 0.25 nm give at T = 1400 K about a quarter 
of the msqd. At T = 1000 K this contribution is dimin- 
ished to 10%. This means that at T = T c short jumps 
dominate the diffusion of both components. 




time (ns) 



FIG. 10: Average atomic mean square displacement versus 
time at T = 900 K (full line). Atomic mean square displace- 
ment calculated from the jumps used for the distribution, 
Fig. 1, for different cutoffs: from top to bottom all jumps with 
jump lengths greater than 0.06, 0.1, 0.15, 0.2 and 0.25 nm. 

At T = 900 K, in the glass, the situation is similar 
for Cu. For Zr the contribution of long jumps is nearly 
negligible. The contribution of short jumps, i < 0.1 nm 
and particularly £ < 0.06 nm (dashed and dotted lines in 
Fig. nj]J, is severely overestimated. This can be traced 
to the "elastic displacements" discussed in the previous 
section. The motion of the Cu-atoms in the Zr matrix 
causes displacements of the latter. Since the Zr subsys- 
tem is sluggish compared with that of Cu, it has some 
memory over several jumps of the Cu system. Zr-atoms 
will have a preference to return to their original site. Or 
speaking in terms of potential energy, the system can 
be approximated for some time as being in a minimum 
with the Cu motion as perturbation. The Zr-atoms re- 
main in their "cages" over many "cage-escapes" of the 
Cu-atoms. While the Zr is in its cage it will occasionally 
perform jumps picked up by the algorithm, followed by 
several short jumps, not picked up. This can be seen by 
monitoring the motion of the single atoms. An indica- 
tion of this effect can be seen in Fig. 1101 where, for Zr, the 
dotted curve initially follows the real msqd (full line) but 
later rises also when the real msqd of Zr stays more or 
less constant while the msqd of Cu increases. This is the 
discussed effect of Zr-atoms temporarily being displaced 
by Cu-jumps and their return in several smaller steps not 
being detected. In this sense the "surplus" of detected 
Zr-jumps corresponds to an in-cage motion. 



We argued above that the rapid increase of P J (T, £) for 
I — > is due to elastic displacements, i.e. displacements 
caused in the matrix when an atom or a group of atoms 
jumps, and that these displacements will not contribute 
strongly to diffusion. As diffusion drops with sinking 
temperature, the ratio of probability of these non diffu- 
sive "jumps" over the one given by Eg. 1151 rises strongly, 
see Fig. |SJ If one does a simple correction of the curves 
for the two smaller cutoffs (€ cu tt = 0.06 and 0.1 nm) by 
this ratio, the excess is removed and the msqd is actually 
underestimated by 30%. Such a simple correction does, 
of course, not distinguish between "in cage" and "out of 
cage jumps" . These "non-diffusive short range jumps" 
exist also in lattices. A jump of an atom into a neighbor- 
ing vacancy site causes displacements of the surround- 
ing atoms which will disappear again when the vacancy 
moves on. In lattices these can easily be measured and 
accounted for by means of lattice geometry and symme- 
try. In an amorphous material this is no longer possible 
and it becomes only a posteriori clear which of the short 
range displacements contribute to diffusion and which are 
non-diffusive. 




time (ns) 



FIG. 11: Average atomic mean square displacement versus 
time at T = 800 K (full line). Atomic mean square displace- 
ment calculated from the jumps used for the distribution, 
Fig. 1, for different cutoffs: from top to bottom all jumps with 
jump lengths greater than 0.06, 0.1, 0.15, 0.2 and 0.25 nm. 



Decreasing the temperature to 800 K the contribu- 
tion of these elastic or "in cage" jumps increases further 
and becomes visible also for the more mobile Cu-atoms, 
Fig. 1111 This can be attributed to the increasing dynamic 
heterogeneity. Like most of the Zr-atoms, an increasing 
number of Cu-atoms becomes immobile on the time scale 
of several jumps. The Cu-atoms, mobile at a given time, 
move in a matrix of Zr- and immobile Cu-atoms which 
are temporarily displaced. 
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VII. CORRELATION FACTOR 

In analogy to Eq. El we can define an "uncorrelated 
msqd" by 



|Ai?uncorr(i,4ut)| 2 = ( £ |AR"(t') 



(18) 



t' <t 

AS"X r 



and a correlation factor 38 



/ corr = lira 



A /? n 

jump 



(t,4 



(t,4ut)l' 



(19) 



Correlation factors have been studied in the past ex- 
tensively for diffusion in lattices 39 40 . If the jumps are 
completely uncorrelated (random walk) one has / CO rr = 1 • 
This would be the case for an isolated interstitial atom, 
e.g. H in Pd, moving between equivalent sites. For dif- 
fusion by a vacancy mechanism, neglecting lattice dis- 
placements, the correlation factor is reduced in simple 
lattices to / CO rr ~ 1 — 2/z where z is the number of near- 
est neighbors. The term 2/z accounts for the direct back- 
jumps of the tracer atom into the same vacancy. Exact 
calculations give for the vacancy mechanism in fee lat- 
tices / corr = 0.78, whereas for a diamond lattice one has 
/coir = 0.5 due to the lower number of neighbors. In 
simple lattices the correlation factor is given by geome- 
try alone and is temperature independent. In more com- 
plicated structures when more than one type of jump is 
involved the correlation factor becomes temperature de- 
pendent. An bias due to an external field increases the 
correlation factor. 

If one would study self diffusion in a simple fee lattice 
along the lines of the two preceding sections one would 
obtain two sets of jumps with two sets of correlation fac- 
tors. "Diffusive jumps" would have a jump length of 
about a nearest neighbor distance and a correlation fac- 
tor as discussed above. The short "non-diffusive jumps" , 
i. e. the temporary displacements due to the "diffu- 
sive jumps" , on the other hand, would have a correla- 
tion factor / corr = 0. The total correlation factor is thus 
somewhat reduced from the rigid lattice value. Total cor- 
relation factors zero are observed in lattices when atoms 
jump between a finite number of sites, e. g. the cage 
motion of interstitial Co in Ai^i. 

These general considerations assume that the diffus- 
ing particle is completely thermalized between jumps, i. 
e. that there are no memory effects. At temperatures 
near the melting point this is no longer the case. In 
simulations of vacancy diffusion in Al and Na lattices it 
was found, that when the waiting time between succes- 
sive jumps shrinks to the order of vibrational time scales 
(ps), successive jumps become correlated, forward jumps 
become more frequent and, subsequently, the correlation 
factor for vacancy diffusion becomes larger than the ran- 
dom walk value, / corr > l 40 i 42 . 

Transferring these results to the amorphous and under- 
cooled liquid states, one expects to find total correlation 



factors considerably smaller than 1 for diffusion via a va- 
cancy type mechanism, and near to 1 for diffusion via an 
inherent mechanism as postulated in the review 1 -. The 
correlation factor for the longer jumps should be larger 
than for the shorter ones since the latter have a larger 
non-diffusive contribution. Cooling below T g , the corre- 
lation factor should drop because more and more jump 
directions become blocked. Unfortunately, due to the 
computational limitations, this effect cannot be explored 
fully at present. In the liquid, where we have seen that 
the jump frequency becomes comparable to the vibration 
time, we expect an increase of the correlation factor with 
temperature. 
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FIG. 12: Average atomic mean square displacement versus 
time at T = 800 K (full line). Atomic mean square displace- 
ment calculated from the jumps used for the distribution, 
Fig. 1, for different cutoffs: from top to bottom all jumps with 
jump lengths greater than 0.25, 0.20, 0.15, 0.1 and 0.06 nm. 



The correlation factors calculated from the observed 
jumps are shown in Fig. 1121 By their very construc- 
tion immediately following jumps are excluded due to the 
dead time between jumps. The general trends are clearly 
seen in the Fig. The correlation factors for the long jumps 
are close to 1 or even larger, indicating diffusion by an 
inherent mechanism. The total correlation factor drops 
upon cooling which explains the drop of the diffusion 
constant below its Arrhenius line despite the seemingly 
constant effective activation energy for the jumps them- 
selves. Striking is the increase of the correlation factor 
for Cu jumps to values around two for T = 1400. This 
will be an effect of insufficient thermalisation between 
jumps. It also indicates the transition from jump to flow 
motion. For the larger and heavier Zr-atoms this effect 
is strongly suppressed. 
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VIII. DISCUSSION 

Upon quenching, the mean square displacements, in 
a log-log representation, show for both components of 
Cu33Zr67 the normal behavior of a transition from a hot 
liquid, via an undercooled one, to a glass. Between the 
oc t 2 behavior, indicating ballistic motion or vibrations, 
and the oc t regime of diffusion a plateau develops. The 
diffusion dynamics can be described by the mode cou- 
pling scenario. However, the temperature dependences 
of the two components do not run as parallel as in the of- 
ten used binary Lennard- Jones system either under con- 
stant volume^ or constant pressure conditions^. Such 
a difference could be the origin of the order of magni- 
tude difference between diffusivity of Ni and Ti and the 
one gained via the Einstein relation from the viscosity in 
ZrTiCuNiBe bulk glasses^ which was, similarly to the 
present findings, explained by a faster diffusion of the 
smaller components in a relatively rigid Zr-matrix. This 
behavior seems to be typical for ZrCu systems. In a sim- 
ulation of the binary NiZr system no such pronounced 
effect was foundi 

Upon quenching, the dynamic heterogeneity rises 
rapidly above its vibrational value at temperatures be- 
low T = 1400 K which, in absence of a reliable estimate 
of the melting point of our system, we take as a rough 
indication of undercooling. 

The distribution of the jump lengths of both compo- 
nents, both above and below T c is a smooth function 
of distance, definitely excluding a preference of nearest 
neighbor jumps. Such a smooth distribution is to be ex- 
pected for collective jumps as they are seen in the typical 
chain-like motion of metallic glasses&2*ii. Using the en- 
ergy landscape approach similar results were found for 
a binary Lennard-Jones system near and were also 
attributed to cooperative motion which in that system 
would be additionally furthered by the inherent higher 
density. It is encouraging that the two different and com- 
plimentary approaches lead to the same conclusion. 

This result, of a missing typical jump length, will 
not be affected by the rapid quench inherent to simula- 
tions. Both experiment*! and simulation of either activa- 
tion volumeiS or isotope effect indicate that insufficient 
aging would, if anything, enhance the number of single 
particle jumps with their inherent length scale. 

For the structure forming majority component Zr, the 
jump length distribution is given by the product of two 
terms, an Arrhenius term for the temperature depen- 
dence and a temperature independent distribution of 
jump lengths. An Arrhenius law for the temperature 
dependence, irrespective of the details of the barrier dis- 
tribution has been derived earlier in an effective medium 
treatment of hopping in disordered materials^. The spa- 
tial dependence shows a power law for short distances and 
and exponential one for the longer distances. This form 
is again in agreement with the findings of the energy 
landscape study 4 ^. The power law part for short dis- 
tances can be understood as being caused by the elastic 



displacement accompanying any hopping process. These 
would be also seen for hopping in a lattice where, differ- 
ent from a disordered system, they can easily identified 
from translational symmetry. They are largely reverted 
by subsequent hops and will not contribute markedly to 
long range diffusion. The exponential long distance de- 
pendence results, in our opinion, from the collectivity 
of the jumps. Jumps are closely correlated to low fre- 
quency quasi-localized vibrations^. The cores of these 
show, similarly to the true localized vibrations, an expo- 
nential decay of the amplitudes^. 

The jump length distribution of the minority compo- 
nent, Cu, has qualitatively the same behavior, however, 
with one important difference. The Arrhenius scaling 
with temperature breaks down above T c . The distribu- 
tion is stretched to greater length and also the total num- 
ber of long jumps is increased. We have related this to 
lacking thermalisation between jumps, a view supported 
by the correlation factors. 

The correlation factors are of order unity. This is a 
signature of diffusion by an inherent mechanism, as op- 
posed to a defect (quasi-vacancy) mediated one which 
would lead to clearly smaller correlation factors. The 
correlation factor of Cu for T > T c is an exception. It 
increases strongly with temperature. Such correlation 
factors > 1 are known to result when the time between 
jumps no longer allows for thermalisation. It indicates, 
in this respect, the transition from jump-like motion to 
flow. 

From the contribution of the different jump length to 
the mean square displacements it is obvious that long dis- 
tance jumps only play a minor role in Zr diffusion. Their 
role is greater for Cu where at T c jumps over more than 
0.2 nm contribute more than 30%. From direct inspec- 
tion of individual jumps one sees, that in the usual chain 
of collectively jumping atoms Cu-atoms will show the 
largest displacements. At the higher temperatures, often 
more than one Cu atom jumps over a nearest neighbor 
distance. The replacement chains observed in the early 
work 5 ^ are a special case of this scenario. 

The absence of a preferred jump length might, prima 
facie, be taken as a contradiction to the time evolution of 
the van-Hove self correlation function where, in the un- 
dercooled liquid, clearly a secondary peak at the nearest 
neighbor distance evolves. One of the authors showed 
earlier for the same material that the time evolution of 
the van-Hove self correlation function for Cu near T = T c 
can be reproduced by a model comprising only nearest 
neighbor jumps and a residue of short distance jumps, 
called flow-like motion^. For Cu an apparent activa- 
tion energy of 850 meV, in good agreement with the 
present value of 820 meV, was given. Combining this 
and the present results one cane conclude that the sec- 
ondary peak in the van-Hove function indicates preferred 
resting positions for the Cu-atoms in the Zr-matrix which 
is rigid on the time scale of the Cu diffusion. A similar 
conclusion was arrived at for the Lennard-Jones system 
from the energy landscape pictur e 45 ! 50 when the dynamic 
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heterogeneity creates a sufficiently rigid matrix of im- 
mobile particles. Diffusion becomes dominated by tran- 
sitions between these minima, a connected network of 
minima is formed^iiSiS or "coarse graining" occurs^. 
The strongly different diffusivities of the present sys- 
tem enhance the effect and are an additional source of 
"metabasins" 55 . The diffusion of the Cu in the deeply 
undercooled liquid and even more in the glass has fea- 
tures of a diffusion between "traps" . 

In the previous worki&, as the temperature is low- 
ered, an increase in jump reversal was reported. This 
is compatible with the reported correlation factors of the 
present work. A relative increase of reversible jumps for 
lower temperatures has also been observed from a jump 
analysis of the binary Lennard-Jones system below T<^. 



IX. SUMMARY 

In summary we calculated the diffusion of Cu and Zr 
above and below the glass transition temperature. Dif- 
ferent to the commonly studied Lennard-Jones system 
Cu shows in the undercooled liquid a strongly enhanced 
diffusivity, compared to the majority component Zr. We 
think that this effect is related to the weaker coupling be- 
tween the two components in the present system which 
is seen also in the reduced enthalpy of fusion. 

The distribution of atomic jumps in the liquid and 
glassy state can be described by simple exponential de- 
pendencies on jump length and temperature. There is a 
smooth transition from flow to hopping. For the faster 
Cu this is reflected in the undercooled liquid by corre- 
lation factors exceeding unity, indicating a breakdown 
of thermalisation between jumps. We find no preferred 
jump length around the nearest neighbor distance. The 
observed secondary peaks in the van-Hove function are, 
therefore, clearly not an effect of the jumps directly but of 
increased waiting times at given sites. This effect will be 
enhanced by the growth of the dynamical heterogeneity 
upon cooling. On the time scale of the inverse of the jump 
frequency of the mobile atoms, more and more atoms are 
immobile and provide a semi rigid background thus cre- 
ating preferred sites. The effect is seen more strongly by 
the smaller atoms (Cu). 



APPENDIX A: MODIFIED EMBEDDED ATOM 
INTERACTION (MEAM) 

We use the MEAM model developed for the CuZr sys- 
tem by GaukelSI. For completeness we give below the 
analytic form and the parameters. 

In the Embedded Atom Method (EAM) model^ the 
interatomic interaction is described by pair potential 
<I>two and an embedding energy F c& (p) which accounts 
for the additional many body effects due to the electronic 



density: 

N N 

£ P °t= e 2<o"(^ r n+E F c;(p(R")) (ai) 

m,n=l n =l 

where R mn = |R m — R n | is the distance between atoms 
m and n, and l n indicates the type of atom, in our case 
Zr or Cu. 

For the pair potential we use the analytic form 

$ two (ff" n ) = Cl e- C2R " m +c 3 e- c ^ R " m ~ Ro ^ 

+ c 5 (R mn ) 6 +c 6 {K mn ) 7 + c 7 (A2) 

The small parameters C5 , cq and c-j are given by the con- 
dition that "^two and its first two derivatives vanish at 
-Rcutoff- The parameters are given in Tabled for Zr and 
Cu. For the mixed interaction we use the the mean 

4>two U (R mn ) = ($Zo (R mn ) + *£f (R mn )) /2. (A3) 





Zr 


Cu 


ci [eV] 


3.457155 • 10 4 


6.547955 • 10 5 


C2 


4.479563 


6.487234 


c 3 [eV] 


-1.062312 


-2.383022 




0.8614267 


0.1318588 


Ro [A] 


2.952510 


0.6049550 


C5 


eV/A 6 




-4.971823 • 10 -8 


-3.085439 • 10~ 4 


C6 


eV/A 7 




6.337047 • 10~ 9 


5.310256 ■ 10~ 5 


c 7 [eV] 


6.636436 ■ 10~ 4 


0.9001767 


-Rcutoff [A] 


6.57620462 


4.44761582 



TABLE I: The constants of the pair potential, Eq. IA2I 



For the embedding term, we use the form proposed by 
Basket for both components 

Fea(p(R n )) = c lP (R n ) ■ ln(c 2P (R")) (A4) 

with the parameters given in Table ITT1 





Zr 


Cu 


Cl 


eV-A 3 




2.71082977 


2.02891465 


C2 


A 3 


0.737740301 


0.876399217 



TABLE II: The constants used in the embedding term, 
Ea.lAH 
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In the original EAM, the density p in Eq. IA1I and 
Eq. IA4I is given as a superposition of radial functions 



N 



Po 



(R n ) = f( Rmn )- 



(A5) 



For the function f(R mn ) we use an exponential plus ad- 
ditional terms to set the function and its first two deriva- 
tives to zero at the cutoff. 

f(R mn ) = Cl e" C2iri " + c 3 (R rnn ) Ci + c 5 (i? m «) C4+1 + C6 

(A6) 

The parameters are compiled in Table. II I II 





Zr 


Cu 


Cl 


l/A 3 




1.04537268 


0.896989894 


C2 


0.123135389 


0.286315587 


<=3 


l/A 3 




5.54697310 ■ 10~ 12 


6.74032288 ■ 10~ 6 


C 4 


15.1751328 


7.51090487 


C5 


l/A 3 




-9.43929033 • 10" 13 


-1.22957801 • 10" 6 


C6 


l/A 3 




-0.597896189 


-0.345100552 


-Rcutofi A 


5.16447906 


4.28405051 



TABLE III: Parameters of the spherical density function, 
Ea.lXel 



In the MEAM additional angular terms are added to 
allow for covalent effects^ needed to describe non ideal 
hep lattices. We restrict ourselves to terms in the third 
power of the cosine of the apex angle 



cos(9 m7li ) 



R mn R h 

The density then takes the form 
p(R n )=Po (R")exp 



(po(R")) 2 

cos 3 (6 mn/ ) ff(R mn ) fi{R ln ) 



C3 



N 



(A7) 



The angular correction in the MEAM is only needed 
for Zr as apex atom but not for Cu. We, therefore, put 
c 3 = for Cu, Table 

For the radial function in the angular correction term 
we used the same form as for the spherical part 

f 3 (R mn ) = Cl e- C2Rmn + c 3 {R mn ) C4 + c 5 {R™) c * +1 + c 6 

(A8) 

with the parameters given in Table Ivl 
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